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1. OVERVIEW 

This is a final report as far as our work at University of Minnesota is concerned. The report 
describes our research progress and accomplishments in development of high performance 
computing methods and tools for 3D finite element computation of aerodynamic characteristics and 
fluid-structure interactions (FSI) arising in airdrop systems, namely ram-air parachutes and round 
parachutes. This class of simulations involves complex geometries, flexible structural components, 
deforming fluid domains, and unsteady flow patterns. The key components of our simulation 
toolkit are a stabilized finite element flow solver, a nonlinear structural dynamics solver, an 
automatic mesh moving scheme, and an interface between the fluid and structural solvers; all of 
these have been developed within a parallel message-passing paradigm. 


2. METHODS 

The methods used in this project include efficient flow solvers based on the space-time finite 
element method and parallel computing strategies using the message-passing model,. In parachute 
systems, where the deformations are large, the structure is currently modeled by cable-membrane- 
lumped mass elements. The nonlinear response is accounted for by casting a Lagrangian 
description of the governing equations. In computation of fluid-structure interactions involved in 
parachute systems, we are increasing the sophistication of handling the coupling between the 
aerodynamics and the parachute deformations. Also, more sophisticated mesh generation and 
update methods are to be developed to handle various stages of the ram-air parachute, such as 
turning with unsymmetric flare. Implementation of these methods is being provided on a wide 
variety of architectures such as the CRAY T3E and SGI MP systems. 


3. FSI SIMULATIONS AND RESULTS 

In fluid-structure interactions, we developed a parallel cable-membrane structural dynamics solver. 
A coupling strategy for structural and flow solvers was also designed; this accommodates both 
nodally equivalent or incompatible meshes. The coupled solver was used to simulate a time- 
accurate fluid-structure interactions for a T-10 personnel parachute, and also a steady-state 
configuration of a ram-air parachute. In the T-10 simulations, the canopy is observed to inflate in 
certain regions and deflate in other regions, concurrent with the pressure fluctuations on the canopy 
surface due to vortex shedding. In the ram-air parachute, the initial configuration of the parachute 
is a smooth curved wing. Upon inflation we observe the appearance of "bumps" and a change in 
its spanwise curvature. 

We tested our methods on time-accurate simulation of flow past round and cross parachutes at 
realistic conditions using unstructured meshes. The geometry for the cross parachute was obtained 
from researchers at Natick RDEC and the computed drag coefficients are in good agreement with 



experimental results. The capability of our automatic mesh generator was enhanced to create dual 
surfaces on the parachute canopies, since these are required when modeling them as structural 
membranes. A parallel cable-membrane structural solver was developed and interfaced with a 
space-time based incompressible flow solver. Preliminary results for round and ram-air parachutes 
were omained using the coupled solvers. In addition, stand alone structural dynamics simulations 
for the X-38 ram-air parafoil were conducted at the AHPCRC and Natick RDEC (with assistance 
from the University of Connecticut). These models are utilizing the "cut pattern" geometry of the 
parafoil, internal ’ribs" and suspension lines. 
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1. INTRODUCTION 

We are interested in developing a fluid-structure interaction (FSI) capability 
which can be used for 3D simulation of airdrop systems. The systems under com 
sideration include conventional personnel round parachutes, cross parachutes which 

can lmitedghde f Capab i 1 , ty ? nd large gliding ram ’ air parachutes (parafoils) which 
can carry upto 21 tons. The design emphasis for these systems is precision dehverv 

hen deployed under demanding conditions such as strong wind gusts and large 
offsets. Since conventional design techniques are time consuming Expensive and 

Et ' 0UI aim iS t0 P rovide a ™ble, cost effective design aid 
on high performance computing technology. 

"" SyStemS P resent y er y complex dynamics arising from interactions be- 
tween the canopy, suspension lines, payload and surrounding air. Parachutes can 
experience ^significant canopy deformations and changes in orientation, even Turing 
steady-state operations. Modeling these phenomena involves a time- variant domain 

^ phys^ 57 m mUSt be aCCOUnted for in order to correctly represent 

The early models developed by us [1, 2) represented the canopy and payload as 
a system of point masses where the deformation of the canopy was assuL^d to be 
given. In a complementary effort, Benney et al. [3] present^ detaSTtmctural 
anaiysis of parachute systems with assumed pressure distributions. Stein et al [41 
presented an axisymmetnc model where the parachute was represented by cable 
elements, and its response was coupled to the flow field. ^ 

i'°. mbbles a stabilized space-time finite element flow 
solver w!th a cable-membrane finite element structural dynamics solver. Both solvers 

dsJLwT ^ ° ped mterfaced Wlthin a message-passing parallel computing 
paradigm^ The other key components in our FSI simulation toolkit are a coupling 
strategy between fluid and structural solver, and a mesh mover for deforming the 
fluid mesh based on displacements at the fluid-structure interface. g 

2. FLUID DYNAMICS SOLVER 

2.1 Governing Equations 

Letjl{ C IF* be the spatial domain of interest bounded by boundary r{ at 
any instant t . Here the superscript /stands for the fluid component and is the 


number of spatial dimensions. The Navier-Stokes equations governing incompressible 
flows axe: 


P / (-^- + u-Vu-f)-V-<r / = 0 on Sl{ , ( 1 ) 

V • u = 0 on f l{. (2) 

Here p f , u, f and a f are the density, velocity, body force and the stress tensor, 
respectively. The stress tensor is written as the sum of its isotropic and deviatoric 
parts, and the fluid is assumed to be Newtonian. The boundary T{ is composed of 
(T( )g and (T t )h corresponding to Dirichlet- and Neumann-type boundary conditions 
respectively. A divergence-free velocity field is specified as initial condition. A 
Baldwin- Lomax turbulence model is used. 

2.2 Finite Element Formulation 

The flow solver is based on the Deforming-Spatial-Domain/Stabilized-Space- 
Time (DSD/SST) formulation [5, 6] which automatically accounts for deformations 
in the fluid mesh. The formulation is written over individual space-time slabs, and 
successive time slabs are coupled in a weak sense. This allows obtaining the solution 
on a single space-time slab at every instant. The formulation is also stabilized to 
circumvent spurious oscillations due to presence of sharp boundary layers at high 
Reynolds number [7] and choice of velocity and pressure interpolants which otherwise 
violate the Babuska-Brezzi [8] condition. 

3. STRUCTURAL DYNAMICS SOLVER 
3.1 Governing Equations 

Let Q 3 t c IT 1 "' be the spatial domain of interest bounded by boundary P at 
any instant “t”. Here the superscript s stands for the structural component. The 
governing equations for the structure is obtained from the conservation of linear 
momentum: 


= 0 on fi*. (3) 

Here x s is the structural displacement. As is the case for the fluid, the boundary P 
is composed of (T 4 ) a and (r t 4 ) h . 4 

The structure undergoes large deformations leading to geometric no nlinear ities 
The resulting strains are assumed to be small and hence a materially linear elastic 
model is used. To account for the kinematic nonlinearities, the constitutive equations 
are written in the original undeformed configuration, in terms of the 2nd Piola- 
Kirchoff stress tensor S and the Green- Lagrange strain tensor E. The structure is 
composed of cable and membrane elements. The cable elements are assum ed to he 
in a state of uniaxial stress: 


S 11 = E c G n G n E u . 


The membrane elements are assumed to he in a state of planar stress: 
S* = (X m G ij G kl + fjL m [G a G jk + G*G jl ))E kl , 


(4) 

(5) 


where 


\ 2A m /i m 

/ 'm — 


(A m + 2(x m ) 


Here G j are the components of the contravariant metric tensor E is the modulus nf 

S!* ^ “V" and <> are the ume 

to Bathe [9] and PollO] " rang<! (1 ' 2) ‘ F ° r fUrther deta,k ' the reader is referred 

3.2 Finite Element Formulation 

The finite element formulation is derived from the principle of virtual work: 

k p3 lF ■ 6 *’ dQ + io S - 6Ed * = /. t • &*fT + / ni ft • 6x s d£l. (7) 

15 Sur fi fa f, e fc «ftion which also adds to the nonlinearity since it results in 
a follower force field. The left-hand-side terms in Eq. 7 are written in the original 
configuration. Upon discretization using appropriate function spaces, a nonlinear 

system of equations is obtained at each time-step. This is solved using Newton- 
Raphson iterations: 6 


dt 2 dt 


+ KAd* = R, 
C = r?M + CK. 


S the f ma5 ^ matrix, C is an artificial damping matrix which stabilizes the 
structural system K is the stiffness matrix, R is the residual vector, and Ad s is 
the increment m the nodal displacements. The Hilber-Hughes-Taylor fill scheme is 
used for time-marching. L 1 

4. FLUID-STRUCTURE COUPLING STRATEGY 

Let r{* be the common boundary between the fluid and structure domains Cou- 
pling is established by transferring velocity (and displacements) from structure to 
fluid and in return traction from fluid to structure. Hence from the perspective of 
the flow solver, I\ c (T t ) g and for the structural solver r{ 3 C (H)*. In the current 
implementation only the pressure contribution to the interface traction is considered 
iWo types of meshes can be utilized at the interface: 

• Nodally equivalent meshes. There is one-to-one correspondence between fluid 
and structure nodes at the interface. 

• Incompatible meshes. This allows the use of separately designed meshes for 
fluid and structure and is more efficient (since structure meshes are normally 

coarser than fluid meshes). Data exchanges at the interface require least- 
squares projections. 

The two solvers are loosely coupled in the current implementation, and the solu- 
tion strategy proceeds as follows: 


• Begin time step loop 

- Increment time step 

- Form predictors 

- Begin nonlinear loop 

* Transfer updated velocities and displacement to fluid 

* Solve mesh motion 

* Solve fluid component 

* Transfer updated pressures to structure 

* Solve structure component 
— End nonlinear loop 

— Form correctors 

• End time step loop 

5. MESH-MOVING SCHEME 

In the case of the structural domain, the nodal displacements are computed as 
part of the solution process. For the fluid domain, the interior displacements due to 
motion at rf 3 are obtained by using an automatic mesh-moving scheme [12], where 
the fluid mesh is modeled as a psuedo-linear elastic solid in static equilibrium. 

6. PARALLEL IMPLEMENTATION 

Both the fluid and structure solvers are implemented on distributed memory archi- 
tectures using a message-passing paradigm. The parallel methodology within each 
solver is the same, here the meshes are partitioned into contiguous even-sized clus- 
ters of elements, and each cluster is mapped to a processor. Within each solver, a 
two-step communication route is set up between global and element data structures 
for the gather and scatter operations. 

Communication between the solvers occurs through a master processor. The 
surface mesh at the interface is extracted from the fluid mesh and stored on the 
master processor along with the membrane component of the structure mesh (also 
lying on the interface). A communication trace is stored between the interface meshes 
and the associated solver. Data exchange at the interface is accomplished on the 
master processor and subsequently communicated to the other processors within the 
corresponding solver. 

7. NUMERICAL EXAMPLE 

In a preliminary application, we attempt to obtain the equilibrium shape of a ram- 
air parachute placed in a windstream. The structural mesh consists of 4,799 nodes, 
9,454 two-noded cable elements, 4,796 four-noded membrane elements, and results 
in 13,356 equations. The fluid mesh consists of 247,437 nodes and 236,688 elements 
and results in 1,912,876 unknowns for fluid solver and 684,151 unknowns for mesh 
mover. The interface meshes are nodally equivalent. The lines in the parachute 
structure confluence at a fixed pivot point. The two sides of the parachute canopy 
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Figure 1. FSI application: parachute canopy in equilibrium configuration 


are constrained to move in the vertical plane. Since the flow within the canopv is 
not modeled, an internal pressure is imposed to represent the ram-air inflation The 
parafoil ribs are modeled using cable elements and the lines are treated as slender 
cy in ers to compute their contribution to the aerodynamic drag. The initial shape 
of the parafoil is assumed to be a smooth wing with R/b = 0.6 (where R is the radius 
of curvature and b is the span). Upon placement in the windstream and application 
of internal pressure, the parachute inflates to a ‘bumby’ shape and pitches back to 
a configuration where the moment due to line drag is balanced by the aerodynamic 
moment due to the canopy. Figure 1 shows the surface pressure on the parachute 

canopy together with stream tubes in the equilibrium configuration. This simulation 
was carried out on a CRAY T3E. 


8. CONCLUDING REMARKS 

We presented a methodology for parallel computation of 3D FSI problems. The 
methodology allows us to use compatible or incompatible interface meshes, adding 
a significant amount of flexibility. As a preliminary example, the method was ap- 
plied to predict the equilibrium configuration of a ram-air parafoil placed in a wind 
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1. INTRODUCTION 

Parachute systems are deployed from a variety of aircraft in many different en- 
vironments. These systems deploy a deceleration device, usually made of highly 
deformable fabrics, which must decelerate the payload to a survivable velocity be- 
fore ground impact. Fluid-structure interactions (FSI) are involved at all stages 
of airdrop systems performance including at initial deployment, during inflation 
at terminal descent (or gliding/maneuvering for steerable parachutes), and at soft 
landing (i.e., retraction for round parachutes, flared landing for ram-air parachutes) 
The interactions between the parachute system and the surrounding flow field are 
dominant in most parachute operations, and thus the ability to predict parachute 
fluid-structure interactions is of high interest to the Army [1,2]. 

The dynamics of parachutes are complex and difficult to’model. They are gov- 
erned by a highly nonlinear coupling between the structural dynamics of a highly 
deformable parachute system and the turbulent time-dependent flow field surround- 
ing the parachute canopy. Advances in high-performance computing are making 3-D 
flow simulations and coupled fluid-structure computations for parachute problems 
more economical and feasible [3]. 

A strategy for performing 3-D parachute FSI simulations is being developed and 
initial results for application to a round parachute have been obtained. The FSI 
strategy will be broken into three components: the fluid dynamics (FD) solution 
the structural dynamics (SD) solution, and the coupling of the FD and SD along 
the fluid-structure interface. The FD solution utilizes a stabilized space-time finite 
element formulation [4,5] of the time-dependent, 3-D Navier-Stokes equations with 
a zero-equation Smagorinsky turbulence model [6] to compute the flow fields. The 
problem is discretized with an unstructured finite element mesh to allow efficient 
meshing of the arbitrary spatial domains encountered [7]. For the SD solution, the 
equations of motion for the parachute system are solved using a finite element for- 
mulation for a “tension structure” composed of cables and membranes [8]. The 



coupling of the FD with the SD is implemented over the fluid-structure interface, 
which is the parachute canopy surface. For FD and SD meshes with compatible 
sets of nodes defining the parachute surface, the coupling involves communication 
of necessary information between FD and SD surface nodes. For incompatible FD 
and SD meshes, coupling information must be computed with a more sophisticated 
projection algorithm [9], To handle large deformations experienced in a parachute 
FSI, an automatic mesh moving scheme with occasional remeshing of the spatial 
domain is necessary [10] . This coupling approach has already been demonstrated on 
an axisymmetric parachute inflation simulation [11], 


2. TEST PROBLEM: 3-D FSI FOR A ROUND PARACHUTE SYSTEM 

A T-10 personnel parachute system is composed of a 35-foot diameter canopy 
and 30 29.4-foot suspension lines. The lines connect to two confluence points (which 
approximate the connection points for a personnel harness assembly). This system 
is allowed to inflate when the canopy is subjected to a differential pressure of 0.5 
lb/ft . Figure 1 shows the SD model for the initial, nonstressed configuration of 
the parachute system, which has a flat extended skirt canopy with a vent at the 
apex and suspension lines that continue as 30 reinforcements through the parachute 
canopy. The parachute system is represented by linearly elastic materials, which 
have properties and dimensions representative of a T-10. 



Figure 1. SD model in nonstressed configuration. 
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Figure 2. Fully inflated configuration for T-10 model. 

A 3-D tetrahedral mesh was generated for the FD solution using the three-noded 
membrane mesh for the T -10 canopy in its inflated configuration as the surSe m«h 
Canopy surface nodes were multiply defined, with one node for both the Lmner^and 
lower surface. Figure 3 shows the surface.mesh for the canopy and a slice of the 3-D 

Z t b ZlV he Ca "° Py B Initial . U " Steady flow elutions were obtained for flow 
about the fixed canopy configuration at a Reynolds number of 10 7 using a semi- 
discrete formulation of the incompressible Navier-Stokes equations. Figure 4 shows 
the computed velocities and pressures for the flow field about the T-10 cfLpy frozen 
at one instant in time. The FSI simulations are initiated using the stabilized space- 


tun* foinuilation with the fully developed How field for this confirm ut ion 1 1 s< *< I ,ls 
the initial condition. 



Figure 3. T-10 canopy surface mesh and slice of 3-D mesh bisecting canopy. 
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Figure 4. Computed flow field about T-10 model. 
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1. INTRODUCTION 

With the recent advances in flow simulation and modeling methods and the ad- 
vanced parallel supercomputers, we now have powerful computational tools capable 
of simulating complex flow problems in applied fluid mechanics. Methods designed 
for a general class of challenging flow problems sometimes need to be re-designed 
and/or optimized for more specific classes of problems which pose their own specific 
challenges. One of these specific classes of problems is 3D simulation of unsteady 
wake flow generated by a primary (leading) object and how the wake flow affects a 
secondary (trailing) object placed in this wake. Examples of this class of applications 
are flow around a small aircraft in the wake of a larger aircraft, and flow around a 
parachute crossing the wake flow of an aircraft. In many cases, because the two ob- 
jects are separated by a large distance compared to the length scales of the objects 
these wake regions are rather long, and we need to take this fact into account in 

developing computational methods which will be effective for this particular class of 
problems. 

In this paper, we present a multi-domain parallel computational method for simu- 
lation of unsteady flow past the primary object, unsteady wake flow in the long region 
between the two objects, and the influence of this wake flow on the secondary object. 
The base method is a finite element formulation with the streamline-upwind/Petrov- 
Galerkin (SUPG) [1] and pressure-stabilizing/Petrov-Galerkin (PSPG) [2] stabi- 
lizations. With these stabilization techniques, simulations can be carried out for 
flows with high Reynolds numbers and thin boundary layers, without generating nu- 
merical oscillations but also without introducing excessive numerical dissipation to 
the computations. Furthermore, these stabilization techniques allow us to use equal- 
order interpolation functions (such as trilinear-trilinear and linear-linear) for velocity 
and pressure without generating any numerical oscillations that would normally be 
caused by such combinations of interpolation functions. 

The spatial discretization of the finite element formulation leads to a set of cou- 
pled, nonlinear, ordinary differential equations which are solved with a predictor/mul- 
ti-corrector time marching scheme. Attach time step, a set of coupled, nonlinear 
equations is solved with the Newton-Raphson iterations. At each Newton-Raphson 
r 3, ° f cou P led » linear equations is solved with iterative methods with the 
GMRES [3] update technique. In computation of vector-matrix products involved in 


this innermost iterations, we have been using element-matrix-based, element-vector- 
based (also called matrix-free) [4], and sparse- matrix- based [5] methods. These 
methods, particularly the last two, can be very effectively used for simulation of 
problems with millions of equations. The overall formulation and the solution tech- 
niques have been implemented on parallel computing platforms by using the MPI 
programming environment. The parallel platforms used include CRAY T3D CRAY 
T3E and the SGI POWER CHALLENGE with the R8000 processor. 

The multi-domain computation approach is based on dividing the entire simu- 
lation domain into an ordered sequence of overlapping subdomains. The flow field 
computed over a eading subdomain is used in specifying the inflow boundary con- 
ditions for the following subdomain. The inflow conditions for the first subdomain 
are extracted from the free-stream conditions. Subdomain- 1 is used for computa- 
tion of the unsteady flow around the primary object. Because the primary object 
typically would have a complex geometry, such as an aircraft geometry, the mesh 
constructed oyer Subdomain- 1 would typically be an unstructured one. Therefore 
the computation over Subdomain- 1 would require a general-purpose finite element 
implementation. Subdomain-2 would be used for computation of the unsteady wake 
flow generated by the primary object, and this subdomain would typically be much 
longer compared to Subdomain- 1. Since Subdomain-2 would not involve any ob- 
jects, the mesh constructed over this domain could be a structured mesh, perhaps 
even a uniform one. A special-purpose finite element implementation recognizing the 
uniform nature of the mesh can be optimized to yield much higher computational 
speeds compared to a general-purpose implementation. In fact, the computation over 
Subdomain- 2 can be accomplished by methods other than the finite element method 
such as the spectral method, that might be more desirable for computations over 
very regular geometries. The computation over Subdomain-3, which would contain 
the secondary object, would typically require, similar to Subdomain- 1, an unstruc- 
tured mesh, and consequently a general-purpose finite element implementation. We 
also note that computations over these different subdomains can be performed on 
different computing platforms and almost in parallel, provided a leading subdomain 
is at least one time step ahead of the following subdomain. 

We need to place the inflow boundary of Subdomain-2 sufficiently ahead of the 
outflow boundary of Subdomain- 1, so that the input to Subdomain-2 is sufficiently 
clear of any truncation effects at the downstream boundary of Subdomain-1. Further- 
more, this input needs to be captured early enough before the solution in Subdomain- 
1 enters the coarser regions of the mesh towards the downstream boundary. Similarly 
the inflow boundary of Subdomain-3 needs to be placed sufficiently ahead of the out- 
flow boundary of Subdomain-2, which in turn needs to be sufficiently distant from 
the secondary object. 


To demonstrate this multi-domain approach, we computed flow problems in- 
volving cylinders and wing-shaped objects. In the case of flow past a cylinder our 
purpose was to verify the validity of some of the concepts used and to capture vortex 
shedding m wake regions very far from the cylinder. Due to space limitations, we 
will not be presenting these results in this paper. In the next section, we present the 
computation with three subdomains, where the first and last subdomains include 
wing-shaped objects. 


2 objects CAL EXAMPLE: wing-shaped primary and secondary 

In this 3D simulation, we compute on SGI POWER CHALLENGE and GRAY 
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Figure 1. Arrangement of the wings and the subdomains. 



Figure 2. Unstructured mesh for the trailing wing 



Figure 3. Iso-surfaces for the leading wing (left) and the wake region (right). Left: 
streamwise component of the velocity (corresponding to 97.5% of the free-stream 
t 5 al ) k SpanW1Se ^ streamwise components of the vorticity (corresponding 



Figure 4. Left: vertical component of the velocity at three vertical planes: the right 
wing tip, the center, and the left wing tip. Right: rolling moment for the trailing 
wing, compared to the case when that wing is in a uniform flow field. 
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SUMMARY 

We present a multi-domain computational method 
for simulation of unsteady flow past a primary ob- 
ject, unsteady wake flow in the long region behind the 
P rin >ary object, and the influence of this wake flow 
on a secondary object. The multi-domain computa- 
tion approach is based on dividing the entire simula- 
tion domain into an ordered sequence of overlapping 
subdomains. The flow field computed over a leading 
subdomain is used in specifying the inflow boundary 
conditions for the following subdomain. Over subdo- 
mains which do not involve any complex-shaped ob- 
jects, we use special-purpose implementation of our 
flow solver to carry out the computation more effi- 
ciently. To demonstrate the potential and power of 
this multi-domain approach, we present 3D computa- 
tional results for flow problems involving cylinders and 
wing-shaped objects. 

INTRODUCTION 

Advanced methods we have developed in recent 
years for flow simulation and modeling and the mod- 
ern parallel supercomputers gave us powerful com- 
putational tools capable of simulating complex flow 
problems in many challenging fluid mechanics applica- 
tions. However, methods designed for a general clmm 
of challenging flow problems sometimes need to be re- 
designed and/or optimized to more effectively address 
specific classes of problems which pose their own spe- 
cific challenges. One of such specific classes of prob- 
lems is 3D simulation of unsteady wake flow generated 
by a primary (leading) object and how the wake flow 
affects a secondary (trailing) object placed at a distant 
location in this wake. Examples of this class of appli- 
cations are flow around a small aircraft in the wake of 
a larger aircraft, and flow around a parachute cross- 
ing the wake flow of an aircraft. In many such cases, 
because the two objects are separated by a large dis- 
tance compared to the length scales of the objects, 
these wake regions have to be relatively long. Conse- 
quently, we need to take this requirement into account 
in developing computational methods which will be ef- 
fective for this particular class of problems. 


The multi-domain computational method described 
in this paper was developed for this purpose. With this 
method, we can compute, with two level of parallelism, 
unsteady flow past the primary object, unsteady wake 
flow in the long region between the two objects, and 
the influence of this wake flow on the secondary ob- 
ject. We start with a stabilized finite element for- 
mulation with the stream] ine-upwind/Petrov-Galerkin 
(SUPG) [ 1 ] and pressure-stabilizing/Petrov-Galerkin 
(PSPG) [2) stabilizations. With these stabilization 
techniques, simulations can be carried out for flows 
with high Reynolds numbers and thin boundary lay- 
ers, without generating numerical oscillations but also 
without introducing excessive numerical dissipation to 
the computations. Furthermore, these stabilization 
techniques allow us to use equal-order interpolation 
functions (such as trilinear-trilinear and linear-linear) 
for velocity and pressure without generating any nu- 
merical oscillations that would normally be caused by 
such combinations of interpolation functions. 

The spatial discretization of the finite element for- 
mulation leads to a set of coupled, nonlinear, ordi- 
nary differential equations which are solved with a 
predictor /multi-corrector time marching scheme. At 
each time step, a set of coupled, nonlinear equations 
is solved with the Newton- Raphson iterations. At 
each Newton-Raphson iteration, a set of coupled, lin- 
ear equations is solved with iterative methods with 
the GMRES [3j update technique. In computation 
of vector-matrix products involved in this inne rmo st 
iterations, we have been using element-matrix-based, 
elements vector-based (also called matrix-free) (4), and 
sparse-matrix-based (5) methods. These methods, par- 
ticularly the last two, can be very effectively used for 
simulation of problems with millions of equations. The 
overaU formulation and the solution techniques have 
been implemented on parallel computing platforms by 
using the MPI programming environment. The par- 
allel platforms used include CRAY T3D CRAY T3E 
and the SGI POWER CHALLENGE with the R8000 
processor. 

In the multi-domain computation approach, we di- 
vide the full simulation domain into an ordered se- 
quence of slightly overlapping subdomains. The in- 



flow boundary conditions for a subdomain is extracted 
from tiie the flow field computed over the preceding 
subdomain. The inflow conditions for the first sub- 
domain are extracted from the free-stream conditions. 
Subdomain-1 is used for computation of the unsteady 
flow around the primary object. Because the primary 
object typically would have a complex geometry, such 
as an aircraft geometry, the mesh constructed over 
Subdomain- 1 would typically be an unstructured one. 
Therefore, the computation over Subdomain- 1 would 
require a general-purpose finite element implementa- 
tion. Subdomain-2 would be used for computation 
of the unsteady wake flow generated by the primary 
object, and this subdomain would typically be much 
longer compared to Subdomain-1. Since Subdomain-2 
would not involve any objects, the mesh constructed 
over this domain could be a structured mesh, per- 
haps even a uniform one. A special-purpose finite 
element implementation recognizing the uniform na- 
ture of the mesh can be optimized to yield much 
higher computational speeds compared to a general- 
purpose implementation. In fact, the computation 
over Suhdomain-2 can be accomplished by methods 
other than the finite element method, such as the spec- 
tral method, that might be more desirable for compu- 
tations over very regular geometries. The computation 
over Subdomain-3, which would contain the secondary 
object, would typically require, similar to Subdomain- 
1, an unstructured mesh, and consequently a general- 
purpose finite element implementation. We also note 
that while for each subdomain the computation is car- 
ried out in parallel, a second level of inherent paral- 
lelism can be exploited by recognizing that computa- 
tions over these different subdomains can be performed 
on different computing platforms and almost in paral- 
lel, provided a leading subdomain is at least one time 
step ahead of the following subdomain. 

We need to place the inflow boundary of 
Subdomain-2 sufficiently ahead of the outflow bound- 
ary of Subdomain- 1, so that the input to Subdomain- 
2 is sufficiently clear of any truncation effects at the 
downstream boundary of Subdomain- 1. Furthermore, 
this input needs to be captured early enough before 
the solution in Subdomain- 1 enters the coarser regions 
of the mesh towards the downstream boundary. Sim- 
ilarly, the inflow boundary of Subdomain-3 needs to 
be placed sufficiently ahead of the outflow boundary 
of Subdomain-2, which in turn needs to be sufficiently 
distant from the secondary object. 

To demonstrate this multi-domain approach, we 
computed flow problems Involving cylinders and wing- 
shaped objects. In the case of flow past a cylinder, 
our purpose was to verify the validity of some of the 
concepts used and to capture vortex shedding in wake 
regions very far from the cylinder. In the case of 
flow past wings in tandem, we compute flow past a 
leading wing in Subdomain- 1, in the long- wake region 


(Subdomain-2), and the effect of this woke flow on a 
trailing wing in Subdomain-3. 

NUMERICAL RESULTS 

Cylinder Wake Computation at Re — 300 wit h 
a Highly Refined Mesh " ~ 

We compute the unsteady flow past a circular cylin- 
der to compare the single- and the multi-domain meth- 
ods. The single-domain computation uses the same 
mesh as Kalro and Tczduyar (6|, and the multi-domain 
computation uses a highly refined mesh (see Figure 1). 

Compared to the single-domain model, the number 
of element in the wake domain is ~3.3 times larger 
in the streamwise direction, ~1.8 times larger in the 
crossflow direction, and 2.0 times larger in the span- 
wise direction. 

Figure 2 shows vorticity iso-surface for these com- 
putations. Figure 3 shows the streamwise component 
of vorticity at the centered horizontal plane. 

These figures clearly show that the solution for the 
multi-domain computation catch greater details of flow 
features, and these features are maintained all the way 
to the outflow boundary. 


Cylinder Long- Wake Computation at Re = 140 

The second phase vortex shedding in the far wake 
of a cylinder was first reported by Taneda (7|, followed 
by many other researchers [8-11]. The second phase 
shedding Is reported at He=140, which is within lami- 
nar regime. Cimbala [9] attributes this phenomena to 
hydrodynamic instability, while Williamson [1 1 1 sug- 
gested that this is caused by the sensitivity of the wake 
to small-scale perturbations. 3D numerical simulation 
of this problem is not an easy task because of the long 
domain required (300d). This motivated us to simu- 
late it using a multi-domain computation. 

First a 2D simulation was performed. The entire 
domain was divided into two subdomains: one around 
the cylinder, and the other one in the wake domain, 
starting at 2 i downstream from the center of the cylin- 
der until 300 d, with one element in span wise direction 
(see Figure 4 ). The wake domain model consists of 
1,824,066 nodes and 905,920 hex&hedral elements, re- 
sulting in 5,446,832 nonlinear equations. 

Figure 5 shows the magnitude of vorticity in the 
wake subdomain. No second phase vortex shedding is 
observed In the far wake, though the K ax man vortex 
street is captured in the near wake. From laboratory 
experiments, the transition between these two phases 
is expected at 100~1 50 i. However in our simulation, 
vorticity decays in that region (see Figure 5). This 
leads us to believe that the 3D effects are the cause of 
the second phase vortex shedding. 
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A 3D model is considered as the second computa- 
tion for the same problem with three subdomains: one 
around the cylinder; another in the first wake subdo- 
main, starting at 2d downstream from the center of the 
cylinder until 155d to catch transition; followed by sec- 
ond wake subdomain starting at 150d until 300d. Each 
wake subdomain model has more than four million 
nodes and hexahedral elements, and results in more 
than 17 million nonlinear equations. The computa- 
tion over each wake subdomain requires, for each time 
step, about 36 seconds on 128 processors of the CRAY 
T3E-1200. Figure 6 shows the magnitude of vorticity 
at horizontal center plane in both subdomains. 

The first wake subdomain captures the transition of 
vortex shedding near outflow boundary, and the sec- 
ond one successfully captures the second phase vortex 
shedding. The spacing between the vortices in the 
second phase vortex shedding is about twice as large 
compared to the spacing we see in the first phase. Fig- 
ure 7 shows the streamwise component of the vorticity 
in the first wake subdomain. 

Figure 7 shows the existence of streamwise vortices 
which are similar to the those seen at Re- 300. The 
streamwise vortices interact with the spanwise vortices 
and break down the first phase vortex shedding. From 
these results, we observe that the existence of stream- 
wise vortices, which is a 3D behavior. These vortices 
make the first phase vortex shedding transform to the 
second phase vortex shedding at Re=140. 

Flow Past Wings in Tandem 

In this 3D simulation, we compute on SGI POWER 
CHALLENGE and CRAY T3E unsteady flow past a 
leading large wing and two trailing small wings placed 
in the far wake of the larger one and affected by its 
wing tip vortices (see Figure 8). The Reynolds num- 
ber is 1000 for the leading wing. We assume sym- 
metry with respect to the plane passing through the 
middle of the primary wing and the two trailing ones. 
Subdomain-1 contains half of the primary wing and is 
handled with a general-purpose finite element imple- 
mentation (see Figure 8). Subdomain-2 is the wake 
region and is handled with a special-purpose imple- 
mentation and a structured mesh. Subdomain-3 con- 
tains one of the trailing wings and is handled with a 
general-purpose implementation and an unstructured 
mesh (see Figure 9). The leading and trailing wings 
both have rectangular shapes, an aspect ratio of 8, 
NACA0012 cross-section, and an angle of attack of 8.0 
degrees. The trailing wing has half the cord-length of 
the leading wing. 

Figure 10 shows the results for the leading wing and 
the wake region. The results for the trailing wing are 

‘The mesh generator waa developed by the Team for Ad- 
vanced Flow Simulation and Modeling (T*AFSM) at the Army 
HPC Research Center 


shown in Figure 11. For the leading wing, we see the 
tip vortices as well as a vortex shedding that quickly 
dissipates downstream because of the coarse mesh. For 
the wake region, we see the wing tip vortices and the 
vortex shedding. Figure 11 (left) shows, for the trailing 
wing, the vertical component of the velocity at three 
vertical planes: the right wing tip, the center, and the 
left wing tip. The left wing tip has more exposure to 
the wing tip vortices from the leading wing. Figure 11 
(right) shows the rolling moment for the trailing wing, 
compared to the case when that wing is in a uniform 
flow field. 
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Figure 2. Cylinder wake computation at Re = 300 with a highly-refined mesh. Vorticity iso-surfaces corre- 
spondjng to 0.6 value, obtained with the single-domain computation (left) and the multi-domain computation 
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Figure 5. 2D cylinder long-wake computation at Re = 140. Magnitude of the vorticity for the entire Subdomain- 
2 (upper) and upstream half of Subdomain-2 (lower). 



Figure 6. 3D cylinder long-wake computation at Re = 140. Magnitude of the vorticity for Subdomain-2 
and Subdomain-3 (upper left and right), Subdomain^ (middle), and Subdomain-3 (with smaller magnitude of 
contours) (lower). 



Figure 7. 3D cylinder long- wake computation at Re = 140. Streamwise component of the vorticity at the 
centered horizontal plane, shown for Subdomain-2. 
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Figure 10. Flow past wings in tandem. Left: iso-surfaces for the streamwise component of the velocity for the 
leading wing (corresponding 97.5% of the free-stream value). Right: iso-surfaces for the wake region; streamwise 
and spanwise components of the vorticity (corresponding to 0.5 value). 



Figure 11. Flow past wings in tandem. Left: vertical component of the velocity at three vertical planes: the 
right wing tip, the center, and the left wing tip. Right: rolling moment for the trailing wing, compared to the 
case when that wing is in a uniform flow field. 
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